ANEOS PYROLITE MODEL DOCUMENTATION NOTEBOOK

ANEOS PYROLITE
RELEASE VERSION: SLVTv0.2G1
RELEASE DATE: 20220301

This Jupyter notebook provides documentation for an ANEOS model for pyrolitic melt composition.

Reference: S.T. Stewart, B.A. Chidester, R. Caracas, J. Badro, M.L. Harwell, M. Huff, S.B. Jacobsen, P. Kalita, D.K. Spaulding, A hydrocode EOS for pyrolitic mantles and magma oceans. Lunar and Planetary Science Conference, 53, Abstract 1535, 2022.

When this verion of the model EOS is used, cite this reference and the doi-tagged Zenodo archive that corresponds to the released version of the EOS model on github. The README file in the github repository provides the history of release versions and doi numbers.

Future updates to this EOS will be archived in Zenodo doi-tagged version-controlled releases.

We use two tabulated quality parameters to help guide usage of this EOS model, as described in https://github.com/ststewart/aneos-forsterite-2019/.

Sarah T. Stewart
sts@ucdavis.edu
June 15, 2021
June 20, 2021, v0.1
June 27, 2021, v0.2
March 1, 2022, v0.2 documentation, data comparison updates and official release

VERSION INFORMATION
Forsterite EOS Version: Pyrolite-ANEOS-SLVTv0.2G1
SLVT# translates to solid-liquid-vapor-thermal model adjustment-version number
G1 = the gridded EOS tables used the gridflag=1 definition for the grid points
SLVTv0.2G1 Release Date: 20220301
GitHub: https://github.com/ststewart/aneos-pyrolite-2022/

This is a solid-liquid-vapor ANEOS model with ionization for a simplified pyrolitic silicate composition (Ca$_{0.87}$Fe$_{2.03}$Mg$_{20.22}$Al$_{1.98}$Si$_{16.27}$O$_{58.63}$). This model EOS uses the new heat capacity adjustment factor described in Stewart et al. (AIP Conference Proceedings, 2020) and Stewart, ANEOS Code Modification: Thermal model adjustment parameter, https://github.com/ststewart/aneos-forsterite-2019/EOS-docs/, 2019.

This Jupyter notebook was used to develop and document ANEOS input parameters that are optimized for the liquid and vapor regions. This model EOS is intended for problems involving bodies with bulk silicate Earth composition in the melt and and partial vaporization regimes.

This model was developed with unpublished data (Caracas and Stewart, in prep., and Chidester et al., in prep.), which will be added to the github repository upon publication.

A reference list is provided at the end of this notebook. PDFs of some of the open access technical references are included in https://github.com/ststewart/aneos-forsterite-2019/EOS-docs/.

USAGE
This notebook interacts with multiple outputs from ANEOS, including phase boundary information, the calculated Hugoniot, and tabulated EOS. The notebook also reads in material properties data and extracted data from previous versions of ANEOS model parameter sets for forsterite.

This notebook is designed to be used in two modes:

This notebook is distributed in release mode. Users who wish to use the development mode can obtain the ANEOS source code from: Thompson, S. L., Lauson, H. S., Melosh, H. J., Collins, G. S., & Stewart, S. T. (2019, November 1). M-ANEOS (Version 1.0). Zenodo. http://doi.org/10.5281/zenodo.3525030.

Development mode:
The notebook calls a local aneos executable which expects a formatted ANEOS.INPUT ascii file and an ascii gridtable.txt definition file that is generated by this notebook. ANEOS produces an ANEOS.OUTPUT ascii file and gridded ascii SESAME-format EOS tables.
The user edits (1) the input parameters in ANEOS.INPUT (separately from this notebook) and (2) the output EOS table parameters in the notebook cell below.
This notebook reads some parameters from ANEOS.INPUT and assumes that the variables are spaced by 10 columns.
When the notebook is run, the code issues a system call to run the aneos program in the local directory and reads in results from the ANEOS.OUTPUT, NEW-SESAME-STD.TXT, and NEW-SESAME-EXT.TXT files to make the plots below and to generate the GADGET format EOS table and other miscellaneous files.

Release mode:
This notebook reads in the provided ANEOS.INPUT, EOS tables, and other distribution data to generate the plots below.

OUTPUT/PROVIDED DATA FILES
ANEOS.OUTPUT: Plain text file with an overview of the calculated equation of state. Look at this file for any errors in the EOS calculation.

NEW-SESAME-STD.TXT: Standard length Sandia-style SESAME file with 201 table and 301 table (density, temperature, pressure, sp. internal energy, Helmholtz free energy). 301 table units: g/cm$^3$, K, GPa, MJ/kg, MJ/kg.

NEW-SESAME-EXT.TXT: SESAME-style table with extra variables from ANEOS. Contains the standard 201 table and non-standard 301-extra-variables EOS table. The 301 table has: density grid values, temperature grid values, sp. entropy(T,rho), sound speed(T,rho), sp. heat capacity(T,rho), KPA flag(T,rho). 2-D arrays list all densities, looping over each temperature. 301 table units: g/cm$^3$, K, MJ/K/kg, cm/s, MJ/K/kg, integer flag, integer flag. The KPA flag is an ANEOS output with phase information.

NEW-SESAME-HUG.TXT: Ascii table with the reference state Hugoniot.

NEW-SESAME-MDQ.TXT: Ascii 301-style SESAME table with the MDQ rating.

NEW-SESAME-IEP-BILINEAR.TXT: Standard length Sandia-style SESAME file with 201 table and 301-style table containing the IEP values for bilinear interpolation. Format: 1-D interpolated density points, 1-D interpolated temperature points, 2-D IEP arrays for pressure, sp. internal energy, Helmholtz free energy, sound speed). IEP units are ratios to the original calculated ANEOS point values. Located in the interpolation subdirectory.

NEW-SESAME-STD-NOTENSION.TXT: Standard length Sandia-style SESAME file with 201 table and 301 table (density, temperature, pressure, sp. internal energy, Helmholtz free energy) where the ANEOS tension region is removed and replaced with the solid-vapor coexistence region. 301 table units: g/cm$^3$, K, GPa, MJ/kg, MJ/kg.

NEW-SESAME-GADGETINIT-NOTENSION.TXT: Standard length Sandia-style SESAME file with 201 table and 301 table (density, temperature, pressure, sp. internal energy, specific entropy) where the ANEOS tension region is removed and replaced with the solid-vapor coexistence region. 301 table units: g/cm$^3$, K, GPa, MJ/kg, MJ/K/kg.

NEW-GADGET2-STD-NOTENSION.TXT: Standard GADGET2 EOS table generated by interpolating the ANEOS table with no tension. Format: number of density points, number of sp. entropy points, density grid values, sp. entropy grid values, 2-D arrays of pressure, temperature, specific internal energy, sound speed. 2-D array values list all densities, looping over each value for specific entropy. Units: g/cm$^3$, erg/K/g, dynes/cm$^2$, K, erg/g, cm/s.

NEW-GADGET2-EXT-NOTENSION.TXT: Extended variables GADGET2 EOS table generated by interpolating the ANEOS table with no tension. Format: number of density points, number of sp. entropy points, density grid values, sp. entropy grid values, 2-D arrays of Helmholtz free energy, KPA flag, MDQ flag. 2-D array values list all densities, looping over each value for specific entropy. Units: g/cm$^3$, erg/K/g, erg/g, integer flag, integer flag.

The SESAME tabular equation of state formats are described in Lyon and Johnson (1992).

APPLICATIONS AND LIMITATIONS

Pyrolite is a compositional model for the chemistry and mineralogy of the entire Earth's mantle that approximates the bulk silicate Earth of today, after separation of the iron core (McDonough and Sun, 1995). This composition is based on synthetic glass samples prepared by J. Badro for shock Hugoniot measurements (Chidester et al. AGU Fall Meeting 2021; Stewart et al. LPSC 2022). The initial bulk density of the glass was 2.94 g/cm$^3$. The composition is similar to the one used by Solomatova and Caracas (2019) and Caracas and Stewart (in prep.) for ab initio calculations of the fluid under terrestrial magma ocean conditions and around the vapor curve: NaCa$_2$Fe$_4$Mg$_{30}$Al$_3$Si$_{24}$O$_{89}$. The ab initio data was used as part of this model development.

The reference state is the solid at STP. In ANEOS, the liquid model is a modification of the solid model. As a compromise to span the whole phase diagram, the solid does not represent the true multi-mineral physical properties of a pyrolitic composition. In particular, the solid bulk modulus is lower than the true material to better fit the liquid region experimental data; the increased compressibility partially compensates for the lack of denser, high-pressure minerals. The heat capacity adjustment parameter is fitted to the available shock temperatures and model estimates for the equation of state. The Debye temperature is fitted to attain the desired specific entropy at melting with a single solid phase. Forsterite data was used as a representative of the single component treatment of pyrolitic composition when necessary, e.g. to define the melting point and shock temperatures from an initial solid phase.

This EOS table is appropriate for problems where bulk silicate Earth composition is present as a warm solid near the melt curve, liquid, and partially vaporized. The temperatures in the liquid field are much improved with the adjusted heat capacity compared to previous ANEOS models using the Dulong-Petit limit.

This model does not include a high-pressure phase transition to avoid the artificial discontinuities introduced in the liquid field by ANEOS. This model becomes increasingly erroneous for shock pressures above about 2000 GPa as the curvature on the Hugoniot and the Grueneisen parameter are not correct at extremely high densities (this model does not approach the Thomas-Fermi limit).

This model is not appropriate for problems entirely in the solid mantles of planets. The modulus of the solid is too low and high-pressure polymorphs are neglected.

The functional form for the Gruneisen parameter in ANEOS has difficulty matching the properties of silicate liquids.

The interpolation procedure to make the GADGET2 density-entropy-grid table uses linear interpolation on the same density grid as the density-temperature SESAME table to obtain the EOS values at each GADGET2 entropy grid point. This method introduces substantial errors in the vapor dome.

The calculated critical point pressure is slightly high compared to the DFT estimate and other silicate critical points. The density and temperatures fall in the range constrained by DFT calculations. The vapor curve follows the 50% condensation point for the bulk silicate Earth composition calculated in Lock et al. JGR 2018. Thus this single phase treatment of the multicomponent system provides an estimate for the boiling point of 50% of the mass.

MODEL DEVELOPMENT QUALITY (MDQ) RATING
We provide an MDQ rating which is saved as a 301-style SESAME table. The MDQ rating is a qualitative guide for users that provides some known quality control information. This information is necessarily subjective and incomplete. The rating does not guarantee any fidelity in the EOS model except for points that have been directly compared to experimental data in the plots below. See discussion with the plots below.

INTERPOLATION ISSUES AND THE INTERPOLATION ERROR PARAMETER (IEP)
The IEP measures the quality of the interpolation scheme for a specific variable on a specific density-temperature grid. Robust interpolation of wide-ranging EOS tables is known to be a difficult problem. Most interpolation schemes individually interpolate each variable, which leads to interpolated states that are not thermodynamically self-consistent. While some thermodynamically consistent interpolation schemes have been developed (e.g., Zeman et al. 2019), they more computationally intensive than simpler interpolation methods and have not been widely implemented yet.

The IEP is calculated on the density-temperature grid using the following method:

A tabulated EOS must be sufficiently densely gridded to capture the phase boundaries. Otherwise, non-physical material response can arise during problems that cross the phase boundary. Simple interpolations schemes generally produce large errors in and near phase boundaries. See illustrative plots below.

We recommend that the interpolation scheme, the implementation of the scheme into EOS code functions (e.g., the source code), and IEP values should be reported in published works that use tabular EOS. See the Jupyter Notebook in the subdirectory and the plots below.

ANEOS NOTES
This work uses v1.0 of M-ANEOS released on Zenodo: Thompson, S. L., Lauson, H. S., Melosh, H. J., Collins, G. S., & Stewart, S. T. (2019). M-ANEOS (1.0). Zenodo. https://doi.org/10.5281/zenodo.3525030 This version included Melosh's (2007) treatment for molecular gas and the capability to include a melt curve and solid-solid/liquid-liquid transition (Collins & Melosh LPSC 2014). This version includes an adjustment to the Debye model for the thermal term in the Helmholtz free energy to approach a user-defined heat capacity at high temperatures. The multiplicative factor $f_{cv}$ is entered in input value V44, and the high-temperature heat capacity is $3f_{cv}Nk$.

The ANEOSTEST.f routine was modified to output tabulated EOS. Note that the current version of this function sets positive pressures smaller than 1.E-30 GPa equal to 1.E-30 GPa.
ANEOS2.f was modified to increase the number of points tabulated on the melt curve in the ANEOS.OUTPUT file and to gather the variables for the heat capacity modification.
ANHUG.f was modified to output more Hugoniot points.
ANEOS1.f and ANEOS2.f were modified to increase the high temperature limit for the heat capacity (Stewart et al., AIP Conference Proceedings, 2020).

FUTURE DEVELOPMENTS

CORRECTIONS AND IMPROVEMENTS
Please send corrections to STS and any requests for data to include in the model-data comparison plots.

USER INPUTS FOR SESAME AND GADGET2 TABLE CONSTRUCTION

If the code cell below is hidden, use the button above to reveal the cell.

In development mode, the user must input:

  1. Header information for the SESAME table.
  2. Temperature, density and entropy grid points.

The following code cell also includes the development mode flag and option to skip construction of a GADGET2 table if it is not needed.

Color mesh plots of the SESAME Rho-T table

Model Data Quality (MDQ) Flag

We provide an MDQ rating which is saved as a 301-style SESAME table. The MDQ rating is a qualitative guide for users that provides some known quality control information. This information is necessarily subjective and incomplete. The rating does not guarantee any fidelity in the EOS model except for points that have been directly compared to experimental data in the plots below.

ANEOS Model Compared to Experimental and Ab Initio Data

The following plots compare the ANEOS model Hugoniots with laboratory data and ab initio calculations.

Compared to forsterite Hugoniot (Root et al. 2018) and entropy on the Hugoniot derived by Davies et al., 2019.

Orange points are corrected data points for the liquid (rho=2.597 g/cm$^3$, T$_0$=2273 K) forsterite Hugoniot from Thomas and Asimow 2013.

ANEOS STP Hugoniots: red -- calculated in the ANEOS code for pyrolitic glass density at room temperature interpreted to be in tension; green -- 1 bar glass estimate; blue -- Hugoniot based at STP (model solid at reference density).

Model Critical Point and Shock-Induced Phase Changes

See Caracas and Stewart (in prep.) for the bounds on the critical point for pyrolitic composition.

Model Phase Boundaries

Black lines are the ANEOS phase boundaries. Red curve is the ANEOS Hugoniot. Blue curve is the interpolated Hugoniot from the gridded eos table.

Comparison critical points from experimentally constrained quartz critical point (Kraus et al. 2012), ab initio enstatite critical point (Xiao & Stixrude 2018), ab initio forsterite critical point (Townsend et al. 2020).

Orange diamonds are the 1-bar boiling point; orange squares are the 1-bar melting point for pure forsterite for reference.

Mosenfelder et al. 2007 present an estimate of the melt curve of forsterite at high pressures. For planetary applications, a forsterite EOS is often used to represent the whole mantle. For comparison, the peridotite solidus and liquidus from Fiquet et al. 2010 are shown. The estimates for the mantle solidus and liquidus span a very large range, especially when water is included; see Figure 1 in Andrault et al. 2017.

Temperatures for Di and Di64An36 are calculated in Jing and Karato (2011) Table 10 using P-V shock data from Asimow and Ahrens (2010).

Peridotite composition 3050 K isentrope and liquid density models from Asimow (2018).

ANEOS model comparisons to experimental data

In the liquid field, we assume that the heat capacity is constant and similar to the value for forsterite liquid reported in Thomas and Asimow 2013.

ANEOS Gruneisen parameter and Theta parameter for thermal model

The Gruneisen parameter in forsterite liquid has been determined experimentally from shock wave data (Thomas & Asimow 2013; Davies et al. in prep.). Shockwave data derive the Grueneisen parameter in reference to a known state: \begin{equation} \gamma = \frac{1}{\rho}\frac{(P-P_{\rm ref})}{(E-E_{\rm ref})} \bigg\rvert _V \end{equation}

The ANEOS Gruneisen model is split into two parts. For $\rho \ge \rho_{0}$, \begin{equation} \gamma = \left( \frac{\gamma_0 \rho_0}{\rho} + C_{24} \left( 1-\frac{\rho_0}{\rho} \right)^2 \right) (1-C_{60}) + \left( \gamma_0 + (C_{24}-\gamma_0) \left( 1- \frac{\rho_0}{\rho} \right)^2 \right) C_{60}. \end{equation} For $\rho < \rho_{0}$, \begin{equation} \gamma = C_{16} \rho^2 + C_{17}\rho + 1 + C_{61}. \end{equation}

$C_{16}$ and $C_{17}$ are calculated so that $\gamma$ and $d \gamma / d \rho$ are continuous at $\rho=\rho_0$. The model asymptotes to the Thomas-Fermi limit ($\gamma=2/3$) when $C_{24}=2/3$.

The functional form in ANEOS does not provide a good fit to the data on forsterite. As a result, the model Gruneisen parameter does not extrapolate to the Thomas-Fermi limit at high pressures.

Specfic Heat Capacity and Sounds Speeds

The calculated heat capacities on the forsterite Hugoniot in the liquid region, from Root et al. 2018, are similar to the 1-bar value derived from Thomas and Asimow 2013. Both are larger than the 3nR high temperature value for the solid (Gillet et al. 1991).

Check the grid spacing compared to the phase boundaries

Colors correspond to temperature. Black lines are ANEOS phase boundaries.

Note that ANEOS by default puts a tension region in for the solid that extends to rhomin (V17) in the input file.

If the tension region is kept, it is visible as a kink connects to the low entropy extension of the melt curve in the lower left plot of density vs. specific entropy.

The lack of points below 1000 K in the vapor dome (seen in temperature vs. specific entropy upper left) corresponds to when the sublimation pressure reaches 1.E-30. At this point, ANEOS bottoms out the pressures at this value and fixes the other state variables (it looks like to constant values; will investigate this region more closely later).

The ANEOS triple point convergence is not perfect, and there are kinks in the tabulated phase boundaries right at the triple point.

ANEOS KPA FLAG

                            TABLE          ANEOS
 KPAQQ=STATE INDICATOR      =1, 1p    =1, 1p    (eos without melt)
                            =2, 2p lv =2, 2p liquid/solid plus vapor
                                      =4, 1p solid  (eos with melt)
                                      =5, 2p melt   (eos with melt)
                                      =6, 1p liquid (eos with melt)
                            =-1 bad value of temperature
                            =-2 bad value of density
                            =-3 bad value of material number

To resolve a melt curve, there should be sufficient numbers of points within the melt curve (KPA=5, red points) so that the melt curve is continuously sampled in the grid.

The pink region is where the vapor pressure is bottomed out at 1.E-30 GPa or the pressures are negative (if there is a tension region.)

Every 10th isotherm is plotted.

Check the Melt Curve in Tabulated EOS

Colored lines correspond to isotherms. Black lines are ANEOS phase boundaries.

Tabulations of the equation of state can lead to erroneous material descriptions across the melt curve. In the following plots, the isotherms should be flat (constant pressure) across the phase boundary (black lines). Every 10th isotherm in the table is shown in the plotted regions.

PLOT GADGET2 FORMAT TABLE: Entropy-Density grid

Color contour plots of the GADGET Rho-S table

The SESAME table and GADGET table use the same density array. So the GADGET table is linearly interpolated along a density-column to find the new table values at each point in the GADGET entropy grid.

Errors in the interpolation, particularly near phase boundaries, can be seen in the sound speeds and KPA flags.

References

Andrault, D., Bolfan-Casanova, N., Bouhifd, M. A., Boujibar, A., Garbarino, G., Manthilake, G., ... & Pesce, G. (2017). Toward a coherent model for the melting behavior of the deep Earth’s mantle. Physics of the Earth and Planetary Interiors, 265, 67-81.

Asimow, P. D. (2018). Melts under extreme conditions from shock experiments, in Magmas Under Pressure, Elsevier, p. 387-418. https://doi.org/10.1016/B978-0-12-811301-1.00015-0

Asimow P. D. and Ahrens T. J. (2010) Shock compression of liquid silicates to 125 GPa: the anorthite–diopside join. J. Geophys. Res. 115, B10209. doi:10.1029/2009JB007145.

Asimow P. D., Ghiorso M. S. (1998) Algorithmic Modifications Extending MELTS to Calculate Subsolidus Phase Relations. American Mineralogist 83, 1127-1131.

Caracas, R., Stewart, S. T. (in review) No magma ocean surface after giant impacts between rocky planets.

Centolanzi, F. J., & Chapman, D. R. (1966). Vapor pressure of tektite glass and its bearing on tektite trajectories determined from aerodynamic analysis. Journal of Geophysical Research, 71(6), 1735-1749.

Chase Jr, M. W., Curnutt, J. L., Downey Jr, J. R., McDonald, R. A., Syverud, A. N., & Valenzuela, E. A. (1982). JANAF thermochemical tables, 1982 supplement. Journal of Physical and Chemical Reference Data, 11(3), 695-940. https://janaf.nist.gov

Chidester, B.A., M. Harwell, J. Badro, D. Spaulding, M. Huff, P. Kalita, R. Caracas, S.B. Jacobsen, S.T.Stewart, Shock response of a magma ocean analog to 1 TPa, AGU Fall Meeting, 2021.

Chidester, B.A., Millot, M., Townsend, J.P., Spaulding, D.K., Davies, E.J., Root, S., Kalita, P., Fratanduono, D.E., Jacobsen, S.B., Stewart, S.T., 2021. The Principal Hugoniot of Iron‐Bearing Olivine to 1465 GPa. Geophys. Res. Lett. 48. https://doi.org/10.1029/2021GL092471

Collins, G. S., & Melosh, H. J. (2014). Improvements to ANEOS for multiple phase transitions. 45th Lunar and Planetary Science Conference, 2664.

Costa, G. C., Jacobson, N. S., & Fegley Jr, B. (2017). Vaporization and thermodynamics of forsterite-rich olivine and some implications for silicate atmospheres of hot rocky exoplanets. Icarus, 289, 42-55.

de Koker, N. P., Stixrude, L., & Karki, B. B. (2008). Thermodynamics, structure, dynamics, and freezing of Mg2SiO4 liquid at high pressure. Geochimica et Cosmochimica Acta, 72(5), 1427-1441.

Davies, E.J., Carter, P.J., Root, S., Kraus, R.G., Spaulding, D.K., Stewart, S.T., Jacobsen, S.B. (2019) Silicate Melting and Vaporization during Rocky Planet Formation, J. Geophys. Res. Planets, submitted.

Duffy, T., Madhusudhan, N., & Lee, K. K. M. (2015). 2.07 Mineralogy of Super-Earth Planets. Treatise on Geophysics, Second Edition, Elsevier, Oxford.Duffy, T., Madhusudhan, N., & Lee, K. K. M. (2015). 2.07 Mineralogy of Super-Earth Planets. Treatise on Geophysics, Second Edition, Elsevier, Oxford.

Fiquet, G., Auzende, A. L., Siebert, J., Corgne, A., Bureau, H., Ozawa, H., & Garbarino, G. (2010). Melting of peridotite to 140 gigapascals. Science, 329(5998), 1516-1518.

Gillet, P., Richet, P., Guyot, F., & Fiquet, G. (1991). High‐temperature thermodynamic properties of forsterite. Journal of Geophysical Research: Solid Earth, 96(B7), 11805-11816.

Ghiorso, Mark S., and Sack, Richard O. (1995) Chemical Mass Transfer in Magmatic Processes. IV. A Revised and Internally Consistent Thermodynamic Model for the Interpolation and Extrapolation of Liquid-Solid Equilibria in Magmatic Systems at Elevated Temperatures and Pressures. Contributions to Mineralogy and Petrology, 119, 197-212.

Hosono, N., Karato, S. I., Makino, J., & Saitoh, T. R. (2019). Terrestrial magma ocean origin of the Moon. Nature Geoscience, 1, 418-423.

Jacobsen, S.D., Holl, C.M., Adams, K.A., Fischer, R.A., Martin, E.S., Bina, C.R., Lin, J.F., Prakapenka, V.B., Kubo, A. and Dera, P., 2008. Compression of single-crystal magnesium oxide to 118 GPa and a ruby pressure gauge for helium pressure media. American Mineralogist, 93(11-12), pp.1823-1828.

Jing, Z, and Karato, S. (2011) A new approach to the equation of state of silicate melts: An application of the theory of hard sphere mixtures, Geochimica et Cosmochimica Acta 75, 6780–6802.

Kerley, G. I. (1977). Rational Function Method of Interpolation. Los Alamos National Laboratory, Report LA-6903-MS.

Kraus, R.G., Stewart, S.T., Swift, D.C., Bolme, C.A., Smith, R.F., Hamel, S., Hammel, B.D., Spaulding, D.K., Hicks, D.G., Eggert, J.H. and Collins, G.W. (2012). Shock vaporization of silica and the thermodynamics of planetary impact events. Journal of Geophysical Research: Planets, 117(E9), E09009, doi:10.1029/2012JE004082.

Kumazawa, M., & Anderson, O. L. (1969). Elastic moduli, pressure derivatives, and temperature derivatives of single‐crystal olivine and single‐crystal forsterite. Journal of Geophysical Research, 74(25), 5961-5972.

Lock, S. J. S. J., Stewart, S. T. S. T., Petaev, M. I. M. I., Leinhardt, Z., Mace, M. T. M. T., Jacobsen, S. B. S. B., & Cuk, M. (2018). The Origin of the Moon Within a Terrestrial Synestia. Journal of Geophysical Research: Planets, 123(4), 910–951. https://doi.org/10.1002/2017JE005333

Luo, S. N., Akins, J. A., Ahrens, T. J., & Asimow, P. D. (2004). Shock‐compressed MgSiO$_3$ glass, enstatite, olivine, and quartz: Optical emission, temperatures, and melting. Journal of Geophysical Research: Solid Earth, 109(B5).

Lyon, S. P., & Johnson, J. D. (1992). SESAME: The LANL equation of state database. Los Alamos National Laboratories Report LAUR-92-3407, Los Alamos, NM.

McDonough, W. and Sun, S.-s. (1995). The composition of the Earth. Chemical Geology 120, 223-253.

Marinova, M. M., Aharonson, O., & Asphaug, E. (2011). Geophysical consequences of planetary-scale impacts into a Mars-like planet. Icarus, 211(2), 960-985.

Melosh, H. J. (2007). A hydrocode equation of state for SiO$_2$. Meteoritics & Planetary Science, 42(12), 2079–2098. https://doi.org/10.1111/j.1945-5100.2007.tb01009.x

Mosenfelder, J. L., Asimow, P. D., & Ahrens, T. J. (2007). Thermodynamic properties of Mg2SiO4 liquid at ultra‐high pressures from shock measurements to 200 GPa on forsterite and wadsleyite. Journal of Geophysical Research: Solid Earth, 112(B6).

Mysen, B. O., & Kushiro, I. (1988). Condensation, evaporation, melting, and crystallization in the primitive solar nebula; experimental data in the system MgO-SiO$_2$-H$_2$ to $1.0\times10^{-9}$ bar and 1870 degrees C with variable oxygen fugacity. American Mineralogist, 73(1-2), 1-19.

Nagahara, H., Kushiro, I., & Mysen, B. O. (1994). Evaporation of olivine: Low pressure phase relations of the olivine system and its implication for the origin of chondritic components in the solar nebula. Geochimica et cosmochimica acta, 58(8), 1951-1963.

Richet, P., Leclerc, F., & Benoist, L. (1993). Melting of forsterite and spinel, with implications for the glass transition of Mg$_2$SiO$_4$ liquid. Geophysical Research Letters, 20(16), 1675-1678.

Robie, R. A., Hemingway, B. S., & Takei, H. (1982). Heat capacities and entropies of Mg$_2$SiO$_4$, Mn$_2$SiO$_4$, and Co$_2$SiO$_4$ between 5 and 380 K. American Mineralogist, 67(5-6), 470-482.

Root, S., Townsend, J.P., Davies, E., Lemke, R.W., Bliss, D.E., Fratanduono, D.E., Kraus, R.G., Millot, M., Spaulding, D.K., Shulenburger, L. and Stewart, S.T. (2018). The principal Hugoniot of forsterite to 950 GPa. Geophysical Research Letters, 45(9), 3865-3872.

Solomatova, N. V., & Caracas, R. (2019). Pressure‐induced coordination changes in a pyrolitic silicate melt from ab initio molecular dynamics simulations. Journal of Geophysical Research: Solid Earth, 124, 11232–11250. https://doi.org/10.1029/2019JB018238

Stewart, ANEOS Code Modification: Thermal model adjustment parameter, https://github.com/ststewart/aneos-forsterite-2019/EOS-docs/, 2019.

Stewart, S., Davies, E., Duncan, M., Lock, S., Root, S., Townsend, J., et al. (2020). The shock physics of giant impacts: Key requirements for the equations of state. In AIP Conference Proceedings (Vol. SCCM19, p. 080003). https://doi.org/10.1063/12.0000946

Stewart, S.T., B.A. Chidester, R. Caracas, J. Badro, M.L. Harwell, M. Huff, S.B. Jacobsen, P. Kalita, D.K. Spaulding, A hydrocode EOS for pyrolitic mantles and magma oceans. Lunar and Planetary Science Conference, 53, Abstract 1535, 2022.

Suzuki, A., & Ohanti, E. (2003). Density of peridotite melts at high pressure. Phys. Chem. Minerals 30, 449-456, doi:10.1007/s00269-003-0322-6.

Thomas, C. W., & Asimow, P. D. (2013). Direct shock compression experiments on premolten forsterite and progress toward a consistent high‐pressure equation of state for CaO‐MgO‐Al$_2$O$_3$‐SiO$_2$‐FeO liquids. Journal of Geophysical Research: Solid Earth, 118(11), 5738-5752.

Thompson, S. L., Lauson, H. S., Melosh, H. J., Collins, G. S., & Stewart, S. T. (2019, November 1). M-ANEOS (Version 1.0). Zenodo. http://doi.org/10.5281/zenodo.3525030.

Tillotson, J. H. (1962). Metallic equations of state for hypervelocity impact (No. GA-3216). General Atomics Division, General Dynamics, San Diego, CA.

Townsend, J. P., Shohet, G., & Cochrane, K. R. (2020). Liquid-vapor coexistence and critical point of Mg2SiO4 from ab initio simulations. Geophysical Research Letters, 47, e2020GL089599. https:// doi.org/10.1029/2020GL089599

Vanpeteghem, C. B., Zhao, J., Angel, R. J., Ross, N. L., & Bolfan‐Casanova, N. (2006). Crystal structure and equation of state of MgSiO3 perovskite. Geophysical Research Letters, 33(3), 2005GL024955.

Xiao, B., & Stixrude, L. (2018). Critical vaporization of MgSiO$_3$. Proceedings of the National Academy of Sciences, 115(21), 5371-5376.

Zeman, M., Holec, M., & Váchal, P. (2019). HerEOS: A framework for consistent treatment of the Equation of State in ALE hydrodynamics. Computers & Mathematics with Applications, 78(2), 483-503.

ANEOS references

Collins, Gareth S., and H. Jay Melosh (2014). Improvements to ANEOS for multiple phase transitions. 45th Lunar Planet. Sci. Conf. Abs. 2664.

Melosh, H. J. (2007). A hydrocode equation of state for SiO$_2$. Meteoritics & Planetary Science, 42(12), 2079-2098.

Thompson, S. L. (1990). ANEOS analytic equations of state for shock physics codes input manual. SANDIA REPORT SAND, 89-2951.

Thompson, S. L., & Lauson, H. S. (1974). Improvements in the Chart D radiation-hydrodynamic CODE III: Revised analytic equations of state (No. SC-RR--71-0714). Sandia Labs.

Thompson, S. L., Lauson, H. S., Melosh, H. J., Collins, G. S., & Stewart, S. T. (2019, November 1). M-ANEOS (Version 1.0). Zenodo. http://doi.org/10.5281/zenodo.3525030

Stewart, S. T. (2019). ANEOS Code Modification: Thermal model adjustment parameter. https://github.com/ststewart/aneos-forsterite-2019/EOS-docs/

Stewart, S., Davies, E., Duncan, M., Lock, S., Root, S., Townsend, J., et al. (2020). The shock physics of giant impacts: Key requirements for the equations of state. In AIP Conference Proceedings (Vol. SCCM19, p. 080003). https://doi.org/10.1063/12.0000946

ANEOS Input Parameters

EOS -1: low temperature solid model
THUG, RHUG: initial state for the Hugoniot calculated in ANEOS.OUTPUT.
These values are room temperature and the density of the pyrolitic glass, which ANEOS interprets to correspond to the solid in tension. There is a small difference between this Hugoniot and the Hugoniot for the a glass at room pressure.

V01 nelem = 6: This model has 6 elements: Ca, Fe, Mg, Al, Si, O.
V02 model type=4: solid-liquid-gas model with ionization
V03 rho0=3.35 g/cm3: reference solid density at the reference state, in this case STP.
V04 T0=298 K: reference state is STP.
V05 P0=1.e6 dynes/cm2: reference state STP.
V06 B0=0.95E12 dynes/cm2: bulk modulus is a compromise value to span the phase diagram. The solid is too soft, e.g. the bulk modulus of forsterite at STP is 128 GPa (Kumazawa and Anderson 1969). And the low-pressure liquid is too stiff.
V07 gamma0=0.85: See gamma vs. density plots. This value is OK for the liquid at the reference density.
V08 Tdebye=-1500.0 K: fitted to obtain an entropy of melting similar to pure forsterite. Negative means use full Debye model. True value for forsterite solid is 768K.

V09 TG model=1: based on free-volume theory, generate a derivative of the bulk modulus near 4.
V10 3*C24=4.5: Compromise value for the gamma function to span solid to liquid.
V11 Esep=1.90E11 erg/g: Zero temperature separation energy. Fitted to match the triple point pressure and critical point temperature.
V12 Tmelt=2163 K: Melting temperature at reference pressure. Here, using the value for forsterite when the true system has a solidus and liquidus.
V13 C53=0 erg/g: Critical point adjustment parameter. Not used.
V14 C54=0. [dimless]: Critical point adjustment parameter. Not used.
V15 H0=0: Thermal conduction parameter. Not used.
V16 C41=0: Thermal conduction parameter. Not used.

V17 rhomin=0 g/cm3: minimum density for the solid. Defines the tension region in the model. The default value is 0.8*rho0 when rhomin=0.
V18=0.0 g/cm3: Solid-solid phase transition parameter. Density at onset of transition. Not used.
V19=0.0 g/cm3: Solid-solid phase transition parameter. Density at end of transition. Not used.
V20=0.0 dynes/cm2: Solid-solid phase transition parameter. Pctr, pressure at the center of the transition. Not used.
V21=0.0: Solid-solid phase transition parameter. First derivative of Pctr with respect to density. Not used.
V22=0.0: Solid-solid phase transition parameter. Second derivative of Pctr with respect to density. Not used.
V23 Hfusion=1.45E10 erg/g: Latent heat of fusion at the reference pressure. Fitted to provide a melt curve similar to forsterite.
V24 rhol/rhos=0.90: Volume change on melting. Large volume change similar to 0.9053 for forsterite at 1 bar.

V25 upper=0.0: Upper limit to cold compression curve extension. Default=1. Setting to zero means use default value.
V26 lower=0.: Default=0.
V27 alpha=0.3: Liquid model parameter ($0<\alpha<1$, default=0.3).
V28 beta=0.1: Liquid model parameter ($0<\beta<1$, default=0.1). beta cannot be equal to gamma.
V29 gamma=0.3: Liquid model parameter ($0<\gamma<1$, default=0.2). beta cannot be equal to gamma.
V30 C60=0: Gamma model adjustment parameter. Default=0.
V31 C61=-0.80: Gamma model adjustment parameter to provide best fit to the critical point. Default=0.
V32 C62=0.5: Critical point adjustment parameter. Default=0. ($0<C63<1$). Fitted to improve critical point.

V33 Ionization model=0=Saha model.
V34-V35=0: Reactive chemistry model not used.
Melosh molecular clusters model for the critical point (Melosh fitted values for V36-V43):
V36 Natom=2: number of atoms in molecular clusters
V37 Ebind=4.25 eV: Binding energy
V38 RotDOF=2.0: Rotational degrees of freedom, 2 for diatomic molecule
V39 Rbond=1.5E-8 cm: Length of molecular bond (cm)
V40 VibDOF=1: Numer of vibrational degrees of freedom, 1 for diatomic molecule

V41 Tdebye=2000 K: Vibrational Debye temperature
V42 Mieflag=1: Flag for Mie potential (1) or Morse potential.
V43 aexp=1.7: Power in Mie potential (1 to 2). Fitted to improve critical point
V44 $f{cv}=1.37$: Adjust heat capacity in high temperature limit. $Cv=3f{cv}NkT$.
V45 QCC1: 1.E-30. low density value to transition to ideal gas.
V45 QCC6: 1.E5. high temperature psi value to transition to ideal gas.
V46 to V48=0: Variables not used.

Atomic number and atom fraction for a simplified pyrolitic silicate composition (Ca$_{0.87}$Fe$_{2.03}$Mg$_{20.22}$Al$_{1.98}$Si$_{16.27}$O$_{58.63}$):
8 0.5863
12 0.2022
13 0.0198
14 0.1627
20 0.0087
26 0.0203

Development Notes

STSM 6/15/2021
Set up development notebook for pyrolite composition.
Began with Omega Hugoniot data on pyrolitic glass (PVT) (Chidester et al. in prep.) and DFT-MD calculations by Caracas & Stewart (in review).

STSM 6/20/2021
Initial model (v0.1) is a pyrolitic glass with a forsterite-like melt curve and heat capacity adjustment parameter fitting the shock temperatures.

STSM 6/27/2021
Glass-melt model abandoned because of trouble with bulk modulus drop from solid to melt.
Now, pyrolitic crystalline (dense, STP=3.35 g/cm3) phase with large volume of melting, comparable to 10% volume change on the melting of pure forsterite (from MELTS database). The large volume change leads to a drop in the bulk modulus on melting.
From MELTS: liquid forsterite at 1890°C and 1 bar (Lange): rho = 2687 kg/m^3; solid forsterite at 1890°C and 1 bar (Berman): 2968 kg/m^3. rhol/rhos = 0.90532

Changed compressibility model to TG=1 (free volume theory) which leads to a high-pressure dK/dP of 4, which is what is inferred from the slope of the Us-up relation for pyrolitic glass (in the fluid region).
True Debye temperature for crystalline forsterite is 768+-15 K (Robie et al. 1982). Like for the ANEOS SLVTv1.0 model for forsterite (2019), the Debye temperature value is fitted to produce a nominal specific entropy at the 1 bar melting point (here, chosen to be the forsterite value).
Added comparison plots to data used by Jing and Karato (2011) for their hard-sphere model for silicate melts.
Added comparison plots for mantle model 3050K isentrope and density along the peridotite liquidus from Asimow (2018).

STSM 6/28/2021
Added comparison plots to Misha Petaev's BSE vapor curve from Lock et al. 2018. The pyrolite vapor curve is close to the 50% condensed multi-component curve.
Triple point of forsterite 5.2e-5 bar, 1890 C (2163 K) (Nagahara et al. 1994). Added box for triple phase region from Lock et al. 2018.

More reference data for forsterite as a single phase approximation when needed for the pyrolitic composition: Richet et al. 1993 Enthalpy of melting: T=2174(100) K, dS=464(4.3) J/K/kg. dHmelt=1008736 J/kg = 1.0e10 erg/g
Thomas and Asimow 2013 liquid linear shock Hugoniot: rho0=2.597 g/cm3, T0=2273 K, c=2.67 km/s, s=1.64
Enthalpy of formation at 0 K from JANAF: 1.54E11 erg/g. Referenced to atomic gas at STP.

STSM 3/1/2022
Updating plots and documentation for initial release.

Need to update MDQ and IEP plots.

End of File